#include "BLASTAnalysis.h"
BLASTResult::BLASTResult()
{
}

BLASTResult::BLASTResult(float E, float I, int L)
{
     Thresholds.Evalue_Threshold=E;
     Thresholds.Identity_Threshold=I;
     Thresholds.MappingLength_Threshold=L;
     AmbiCount=0;
     Keywords.push_back("Adlercreutzia");
     Keywords.push_back("Cronobacter");
     Keywords.push_back("Helicobacter");
     Keywords.push_back("Janthinobacterium");
     Keywords.push_back("Rikenella");
     Keywords.push_back("Shuttleworthia");
     Keywords.push_back("unclassified_\"Lactobacillales\"");
}

BLASTResult::~BLASTResult()
{
     //delete[] ReadName;
     //cout<<"1";
     //delete[] TemplateName;
     //cout<<"2";
     DBAbundance.clear(); 
     Abundance.clear();
     //cout<<"3";
     map<string, map<string, MappingInfo> >::iterator it=MappingResult.begin();
     while(it!=MappingResult.end())
     {
         (*it).second.clear();
         it++;
     }
     MappingResult.clear();
     Keywords.clear(); //11-2
     //cout<<"4";
}

int BLASTResult::Import(string fname, GBProcessing& GenBank, RDPIDQuery& RDPQuery, int OperationFlag, int TaxonomicLevel)
{
    string tPrevReadName="N/A",tReadName,tTemplateName;
    map<string,double> LocalAbundance;
    LocalAbundance.clear();
    char dummy[MAX_LINE];
    char** tempList=NULL;
    int HitNum=0;
    float HighScore=0;
    int HighLength=0;
    string BestHit;
	fstream Mapping;
	Mapping.open(fname.c_str());
	if(Mapping.fail()){cout<<"Cannot open BLAST result file"<<endl; return 0;}
    Mapping.getline(dummy, MAX_LINE);
    while(!Mapping.eof())
    {
        if(dummy[0]!='\n'&&dummy[0]!='\0')
        {
            tempList=line_proc(dummy);
            //read in the mapping information
            tReadName.assign(tempList[0]);
            tTemplateName.assign(tempList[1]);
            if(OperationFlag!=0||tReadName!=tPrevReadName)
                read_mapping_info(tReadName, tTemplateName, tempList);
            if(tReadName!=tPrevReadName)
            {
                if(Filter(MappingResult[tReadName][tTemplateName], Thresholds))
                {
                    if(!LocalAbundance.empty())
                    {
                        if(LocalAbundance.size()>1)
                            AmbiCount++;
                        //ambiguity_detect(fname, tPrevReadName, GenBank, RDPQuery, LocalAbundance, TaxonomicLevel);
                        ReadDistr[tReadName]=LocalAbundance; //2012-2-12
                        normalize(LocalAbundance, OperationFlag, BestHit);
                        map<string, double>::iterator itt=LocalAbundance.begin();
                        while(itt!=LocalAbundance.end())
                        {
                            for(int i=0;i<Keywords.size();i++)
                            {
                                if((*itt).first==Keywords[i])
                                {
                                    cout<<tPrevReadName<<"\t"<<Keywords[i]<<"\t"<<(*itt).second<<"\n";
                                    break;
                                }
                            }
                            itt++;
                        }
                        set_Abundance(LocalAbundance);
                    }
                    HitNum=0;
                    LocalAbundance.clear();
                    tPrevReadName=tReadName;
                    BestHits[tReadName]=tTemplateName;
                    BestHit=GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel);
                    HighScore=MappingResult[tReadName][tTemplateName].Score;
                    HighLength=MappingResult[tReadName][tTemplateName].Length;
                    //MappingResult[tReadName][tTemplateName].Score=MappingResult[tReadName][tTemplateName].Identity;
                }
            }
            if(tReadName==tPrevReadName)
            {
                if(OperationFlag!=0||HitNum==0)
                {
                    if(Filter(MappingResult[tReadName][tTemplateName], Thresholds))
                    {
                        if(MappingResult[tReadName][tTemplateName].Score==HighScore)
                        {
                        HitNum++;
                        //MappingResult[tReadName][tTemplateName].Score=(double)(MappingResult[tReadName][tTemplateName].Length/HighLength)*(MappingResult[tReadName][tTemplateName].Score/HighScore)*MappingResult[tReadName][tTemplateName].Identity;
                        if(OperationFlag==2)
                        {
                            map<string, double>::iterator it=LocalAbundance.find(GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel));
                            if(it==LocalAbundance.end())
                                LocalAbundance[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]=MappingResult[tReadName][tTemplateName].Identity;
                            else
                            {
                                if(MappingResult[tReadName][tTemplateName].Identity>LocalAbundance[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)])
                                    LocalAbundance[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]=MappingResult[tReadName][tTemplateName].Identity;
                            }
                        }
                        else
                        {
                            map<string, double>::iterator it=LocalAbundance.find(tTemplateName);
                            if(it==LocalAbundance.end())
                                LocalAbundance[tTemplateName]=1.0;
                            else
                                LocalAbundance[tTemplateName]++;
                        }
                        }//for if HighScore 
                    }//for if(Filter...
                }//for if(OperationFlag!=0||HitNum==0)
            }             
            for(int i=0;i<12;i++)
                DELETE(tempList[i]);
            DELETE(tempList);
        }
        Mapping.getline(dummy, MAX_LINE);        
    }
    if(!LocalAbundance.empty())
    {
        if(LocalAbundance.size()>1)
            AmbiCount++;
        //ambiguity_detect(fname, tReadName, GenBank, RDPQuery, LocalAbundance, TaxonomicLevel);
        ReadDistr[tReadName]=LocalAbundance; //2012-2-12
        normalize(LocalAbundance, OperationFlag, BestHit);
        set_Abundance(LocalAbundance);
    }
    LocalAbundance.clear();
    Mapping.close();
    //MC.greedy_minimum_cover(ReadDistr);
    //MLE.maximum_likelihood_estimate(ReadDistr, MC.SelectedGroups);
    //set_Abundance(MLE.MLEAbund);
    //fstream BestHitList;
    //BestHitList.open((fname+".besthits").c_str(),ios_base::out);
    //map<string,string>::iterator it=BestHits.begin();
    //while(it!=BestHits.end())
    //{
        //BestHitList<<(*it).second<<"\n";
        //it++;
    //}
    //assignment_test(fname, GenBank, RDPQuery);
    cout<<"import "<<fname<<" complete\n";
    return 1;    
}

int BLASTResult::Import_Hierarchy(string fname, GBProcessing& GenBank, RDPIDQuery& RDPQuery, int OperationFlag, int TaxonomicLevel)
{
    string tPrevReadName="N/A",tReadName,tTemplateName;
    map<string,double> LocalAbundance;
    map<string, vector<string> > HierarchyInfo;
    vector<string> tTemplates;
    LocalAbundance.clear();
    char dummy[MAX_LINE];
    char** tempList=NULL;
    int HitNum=0;
    float HighScore=0;
    int HighLength=0;
    string BestHit;
	fstream Mapping;
	Mapping.open(fname.c_str());
	if(Mapping.fail()){cout<<"Cannot open BLAST result file"<<endl; return 0;}
    Mapping.getline(dummy, MAX_LINE);
    while(!Mapping.eof())
    {
        while(tReadName==tPrevReadName&&!Mapping.eof())
        {
            if(dummy[0]!='\n'&&dummy[0]!='\0')
            {
                tempList=line_proc(dummy);
                //read in the mapping information
                tReadName.assign(tempList[0]);
                tTemplateName.assign(tempList[1]);
                read_mapping_info(tReadName, tTemplateName, tempList);
                tTemplates.push_back(tTemplateName);
                for(int i=0;i<12;i++)
                    DELETE(tempList[i]);
                DELETE(tempList);                                       
            }
            Mapping.getline(dummy, MAX_LINE); 
        }
        if(tPrevReadName!="N/A")
        {
        for(int i=1;i<=TaxonomicLevel;i++)
        {
            LocalAbundance.clear();
            HierarchyInfo.clear();
            bool DupFlag=false;
            double Max=0.0;
            int MaxSize=0;
            string MaxIndex;
            for(int j=0;j<tTemplates.size();j++)
            {
                string tTaxon=GenBank.get_Taxonomy(tTemplates[j],i);
                HierarchyInfo[tTaxon].push_back(tTemplates[j]);
                map<string, double>::iterator it=LocalAbundance.find(tTaxon);
                if(it==LocalAbundance.end())
                    LocalAbundance[tTaxon]=MappingResult[tPrevReadName][tTemplates[j]].Score;
                else
                    LocalAbundance[tTaxon]+=MappingResult[tPrevReadName][tTemplates[j]].Score;
                double tSize=LocalAbundance[tTaxon];
                if(tSize>Max){Max=tSize; MaxIndex=tTaxon; MaxSize=HierarchyInfo[tTaxon].size(); DupFlag=false;}
                if(tSize==Max&&tTaxon!=MaxIndex){DupFlag=true;}
            }
            tTemplates.clear();
            for(int j=0;j<MaxSize;j++)
                tTemplates.push_back(HierarchyInfo[MaxIndex][j]);
            if(DupFlag==true)
                cout<<tPrevReadName<<" has >1 best taxonomy at "<<i<<"-th level\n";
            LocalAbundance.clear();
            LocalAbundance[MaxIndex]=1.0;
        }
        }// if tPrevReadName=="N/A"
        if(!LocalAbundance.empty())
        {
            //if(LocalAbundance.size()>1)
            //ambiguity_detect(fname, tReadName, GenBank, RDPQuery, LocalAbundance, TaxonomicLevel);
            normalize(LocalAbundance, OperationFlag, BestHit);
            set_Abundance(LocalAbundance);
        }
        LocalAbundance.clear();
        tPrevReadName=tReadName;
        tTemplates.clear();           
    }
    if(!LocalAbundance.empty())
    {
        //if(LocalAbundance.size()>1)
        //ambiguity_detect(fname, tReadName, GenBank, RDPQuery, LocalAbundance, TaxonomicLevel);
        normalize(LocalAbundance, OperationFlag, BestHit);
        set_Abundance(LocalAbundance);
    }
    LocalAbundance.clear();
    Mapping.close();
    //assignment_test(fname, GenBank, RDPQuery);
    cout<<"import "<<fname<<" complete\n";
    return 1;    
}

int BLASTResult::Import_Pairend(string fname, GBProcessing& GenBank, RDPIDQuery& RDPQuery, int OperationFlag, int TaxonomicLevel)
{
    string tPrevReadName="N/A",tPrevPairName="N/A",tCurrentPairName,tReadName,tTemplateName;
    bool PairFlag=false;
    map<string,double> LocalAbundance1, LocalAbundance2;
    LocalAbundance1.clear();
    LocalAbundance2.clear(); 
    char dummy[MAX_LINE];
    char** tempList=NULL;
    int HitNum=0;
    float HighScore=0;
    int HighLength=0;
    string BestHit;
	fstream Mapping;
	Mapping.open(fname.c_str());
	if(Mapping.fail()){cout<<"Cannot open BLAST result file"<<endl; return 0;}
    Mapping.getline(dummy, MAX_LINE);
    while(!Mapping.eof())
    {
        if(dummy[0]!='\n'&&dummy[0]!='\0')
        {
            tempList=line_proc(dummy);
            //read in the mapping information
            tReadName.assign(tempList[0]);
            size_t found=tReadName.find_last_of("/");
            tCurrentPairName=tReadName.substr(0,found);
            tTemplateName.assign(tempList[1]);
            if(OperationFlag!=0||tReadName!=tPrevReadName)
                read_mapping_info(tReadName, tTemplateName, tempList);
            if(Filter(MappingResult[tReadName][tTemplateName], Thresholds))
            {
                if(tReadName!=tPrevReadName)
                {
                    tPrevReadName=tReadName;
                    BestHits[tReadName]=tTemplateName;
                    BestHit=GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel);
                    HighScore=MappingResult[tReadName][tTemplateName].Score;
                    HighLength=MappingResult[tReadName][tTemplateName].Length;
                    if(tCurrentPairName==tPrevPairName)
                        PairFlag=true;                          
                    if(tCurrentPairName!=tPrevPairName)
                    {                     
                        if(!LocalAbundance2.empty()&&PairFlag==true)
                        {
                            normalize(LocalAbundance2, OperationFlag, BestHit);
                            set_Abundance(LocalAbundance2);
                            PairFlag=false;
                        }
                        HitNum=0;
                        LocalAbundance1.clear();
                        LocalAbundance2.clear();
                        tPrevPairName=tCurrentPairName;
                    }
                }
                if(tReadName==tPrevReadName)
                {
                    if(OperationFlag!=0||HitNum==0)
                    {
                        if(MappingResult[tReadName][tTemplateName].Score<=HighScore)
                        {
                            HitNum++;
                            //MappingResult[tReadName][tTemplateName].Score=(double)(MappingResult[tReadName][tTemplateName].Length/HighLength)*(MappingResult[tReadName][tTemplateName].Score/HighScore)*MappingResult[tReadName][tTemplateName].Identity;
                            if(PairFlag==false)
                            {
                                map<string, double>::iterator it=LocalAbundance1.find(GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel));
                                if(it==LocalAbundance1.end())
                                    LocalAbundance1[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]=1.0;
                                else
                                    LocalAbundance1[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]++;
                            }
                            else if(PairFlag==true)
                            {
                                map<string, double>::iterator it=LocalAbundance1.find(GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel));
                                if(it!=LocalAbundance1.end())
                                {
                                    map<string, double>::iterator it1=LocalAbundance2.find(GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel));
                                    if(it1==LocalAbundance2.end())
                                        LocalAbundance2[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]=LocalAbundance1[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]+1.0;
                                    else
                                        LocalAbundance2[GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)]++;
                                }
                            }
                        }//for if HighScore 
                    }//for if(OperationFlag!=0||HitNum==0)
                }
            }             
            for(int i=0;i<12;i++)
                DELETE(tempList[i]);
            DELETE(tempList);
        }
        Mapping.getline(dummy, MAX_LINE);        
    }
    if(!LocalAbundance2.empty())
    {
        //if(LocalAbundance.size()>1)
            //ambiguity_detect(fname, tReadName, GenBank, RDPQuery, LocalAbundance);
        normalize(LocalAbundance2, OperationFlag, BestHit);
        set_Abundance(LocalAbundance2);
    }
    LocalAbundance1.clear();
    LocalAbundance2.clear();
    Mapping.close();
    //assignment_test(fname, GenBank, RDPQuery);
    cout<<"import "<<fname<<" complete\n";
    return 1;    
}

int BLASTResult::OutputSequence(string oname)
{
    fstream Output;
    Output.open(oname.c_str(),ios_base::out);
    if(Output.fail()){cout<<"cannot open output file\n"; return 0;}
    Output<<Abundance.size()<<"\n";
    map<string,double>::iterator it = Abundance.begin();
    while(it!=Abundance.end())
    {
        Output<<(*it).first<<"\t"<<(*it).second<<"\n";
        it++;
    }
    Output.close();
    return 1;    
}

string BLASTResult::get_SeqID(int pos)
{
    map<string,double>::iterator it = Abundance.begin();
    for(int i=0;i<pos-1;i++)
        it++;
    return (*it).first;
}

int BLASTResult::get_size()
{
    return Abundance.size();
}

double BLASTResult::get_abundance(int pos)
{
    map<string,double>::iterator it = Abundance.begin();
    for(int i=0;i<pos-1;i++)
        it++;
    return (*it).second;
}

bool BLASTResult::Filter(MappingInfo M, Criterion C)
{
     if((abs(M.ReadStart-M.ReadEnd)<C.MappingLength_Threshold)||(M.Evalue>C.Evalue_Threshold)||(M.Identity<C.Identity_Threshold))
         return false;
     else
         return true;
}

int BLASTResult::SetThresholds(float E,float I,int L)
{
     Thresholds.Evalue_Threshold=E;
     Thresholds.Identity_Threshold=I;
     Thresholds.MappingLength_Threshold=L;
}

char** BLASTResult::line_proc(char* dummy)
{
      char** tempList=new char*[12];
      int len=strlen(dummy), count=0;
      len+=(dummy[len-1]=='\n')? -1:0;
      char* temp=new char[len+1];
      memset(temp,0x00,sizeof(char)*(len+1));
      memcpy(temp,dummy, sizeof(char)*len);
      char * pch;
      pch=strtok(temp,"\t");
      while (pch != NULL)                 //read in this line (including mapping destination, similarity, start/end position, etc.
      {
            len=strlen(pch);
            tempList[count]=new char[len+1];
            memset(tempList[count],0x00,sizeof(char)*(len+1));
            memcpy(tempList[count],pch, sizeof(char)*len);
            count++;
            pch = strtok (NULL, "\t");
      }
      delete[] temp;
      return tempList;
}

int BLASTResult::read_mapping_info(string tReadName, string tTemplateName, char** tempList)
{
    MappingResult[tReadName][tTemplateName].Identity=atof(tempList[2]);
    MappingResult[tReadName][tTemplateName].Length=atoi(tempList[3]);
    MappingResult[tReadName][tTemplateName].ReadStart=atoi(tempList[6]);
    MappingResult[tReadName][tTemplateName].ReadEnd=atoi(tempList[7]);
    MappingResult[tReadName][tTemplateName].TemplateStart=atoi(tempList[8]);
    MappingResult[tReadName][tTemplateName].TemplateEnd=atoi(tempList[9]); 
    MappingResult[tReadName][tTemplateName].Score=atof(tempList[11]); 
        /*
        char* pch2;   //extract E-value
        pch2=strtok(tempList[10],"e");
        float tBase=atof(pch2);
        pch2=strtok(NULL,"e");
        float tExponent=atof(pch2);
        MappingResult[tReadName].Evalue=pow(tBase, tExponent);*/
    return 1;
}

int BLASTResult::normalize(map<string,double>& LocalAbundance, int OperationFlag, string BestHit)
{
    double Sum=0.0;
    map<string,double>::iterator it1=LocalAbundance.begin();
    if(OperationFlag!=2)
    {
    while(it1!=LocalAbundance.end())
    {
        (*it1).second=(double)(*it1).second/DBAbundance[(*it1).first];
        Sum+=(*it1).second;
        it1++;
    }
    }//if OperationFlag!=2
    else
    {
    while(it1!=LocalAbundance.end())
    {
        Sum+=(*it1).second;
        it1++;
    }
    } 
    it1=LocalAbundance.begin();
    double Max=0.0;
    double SecMax=0.0;
    string MaxIndex;
    double count=0;
    while(it1!=LocalAbundance.end())
    {
        (*it1).second=(double)(*it1).second/Sum;
        /*
        if((*it1).second==Max)
        {
            //11-10:
            //---------------
            LocalAbundance.clear();
            break;
            //---------------
            //if((*it1).first==BestHit)
                //MaxIndex=(*it1).first;
        }
        */
        if((*it1).second>Max)
        {
            //11-10:
            //----------------
            SecMax=Max;
            //----------------
            Max=(*it1).second;
            MaxIndex=(*it1).first;
        }
        it1++;
    }
    if(OperationFlag==2)
    {
        int a=LocalAbundance.size();
        LocalAbundance.clear();
        //if(a==1)
            LocalAbundance[MaxIndex]=1.0;
    }
    MaxIndex.clear(); 
    return 1;
}

int BLASTResult::set_Abundance(map<string,double>& LocalAbundance)
{
    map<string,double>::iterator it1=LocalAbundance.begin();
    while(it1!=LocalAbundance.end())
    {
        map<string,double>::iterator it2=Abundance.find((*it1).first);
        if(it2==Abundance.end())
        {
            Abundance[(*it1).first]=(double)(*it1).second;
        }
        else
        {
            Abundance[(*it1).first]+=(double)(*it1).second; 
        }
        it1++;
    }
    return 1;
}

void BLASTResult::Iterator_reset()
{
    Iterator=Abundance.begin();
    return;
}

bool BLASTResult::Iterator_end()
{
    if(Iterator==Abundance.end())
        return true;
    else
        return false;
}

void BLASTResult::MappingIterator_reset()
{
    MappingIterator=MappingResult.begin();
    return;
}

bool BLASTResult::MappingIterator_end()
{
    if(MappingIterator==MappingResult.end())
        return true;
    else
        return false;
}

int BLASTResult::MappingIterator_size()
{
    return MappingResult.size();
}

void BLASTResult::DBAbundance_set(string Index, double Key)
{
    DBAbundance[Index]=Key;
    //cout<<DBAbundance[Index]<<"\t";
    return; 
}

int BLASTResult::ambiguity_detect(string fname, string tReadName, GBProcessing& GenBank, RDPIDQuery& RDPQuery, map<string, double>& LocalAbundance, int TaxonomicLevel)
{
    fstream Ambiguity, Ambiguity_B;
    Ambiguity.open((fname+".ambiguity").c_str());
    if(Ambiguity.fail())
        Ambiguity.open((fname+".ambiguity").c_str(), ios_base::out);
    else
    {
        Ambiguity.close();
        Ambiguity.open((fname+".ambiguity").c_str(), ios_base::out|ios_base::app);
    }
    Ambiguity_B.open((fname+".ambiguity_b").c_str());
    if(Ambiguity_B.fail())
        Ambiguity_B.open((fname+".ambiguity_b").c_str(), ios_base::out);
    else
    {
        Ambiguity_B.close();
        Ambiguity_B.open((fname+".ambiguity_b").c_str(), ios_base::out|ios_base::app);
    }
    size_t found=tReadName.find_last_of(".");
    cout<<tReadName.substr(found-8,8)<<"\t"<<RDPQuery.get_RDPID(tReadName.substr(found-8,8))<<"\n";
    string tTemplateName=RDPQuery.get_RDPID(tReadName.substr(found-8,8));
    double Max=0.0;
    string tMaxHit;
    map<string, double>::iterator it=LocalAbundance.begin();
    while(it!=LocalAbundance.end())
    {
        if((*it).second>Max){Max=(*it).second; tMaxHit=(*it).first;}
        it++;
    }
    
    if(GenBank.get_Taxonomy(BestHits[tReadName], TaxonomicLevel)!=GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel))
    {
        Ambiguity_B<<tReadName<<"\t"<<GenBank.get_Taxonomy(BestHits[tReadName], TaxonomicLevel)<<"\t"<<tTemplateName<<"\t"<<GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)<<"\n";                                                                             
    }
    
    if(tMaxHit!=GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel))
    {
    Ambiguity<<tReadName<<"\t"<<tTemplateName<<"\t"<<GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)<<"\n";
    it=LocalAbundance.begin();
    while(it!=LocalAbundance.end())
    {
        Ambiguity<<(*it).first<<"\t"<<(*it).second<<"\n";
        it++;
    }
    }
    Ambiguity.close();
    Ambiguity_B.close();
    return 1;
}

int BLASTResult::assignment_test(string fname, GBProcessing& GenBank, RDPIDQuery& RDPQuery, int TaxonomicLevel)
{
    fstream Assignment;
    Assignment.open((fname+".assignment_test").c_str(), ios_base::out);
    map<string, float> tAssignment;
    map<string, int> tAssignmentWhole;
    vector<string> UniqueMapped;
    map<string, map<string, MappingInfo> >::iterator it=MappingResult.begin();
    int count=0;
    while(it!=MappingResult.end())
    {
        string AssignmentFinal;
        float max=0;
        map<string, MappingInfo>::iterator it2=(*it).second.begin();
        while(it2!=(*it).second.end())
        {
            if((*it2).second.Score>max)
                max=(*it2).second.Score;
            it2++;
        }
        it2=(*it).second.begin();
        while(it2!=(*it).second.end())
        {
            if((*it2).second.Score==max)
            {
                string tTaxon=GenBank.get_Taxonomy((*it2).first, TaxonomicLevel);
                map<string, float>::iterator it3=tAssignment.find(tTaxon);
                if(it3==tAssignment.end())
                    tAssignment[tTaxon]=1;
                else
                    tAssignment[tTaxon]++;
            }
            it2++;
        }
        max=0;
        map<string, float>::iterator it3=tAssignment.begin();
        while(it3!=tAssignment.end())
        {
            if((*it3).second>max)
            {
                max=(*it3).second;
                AssignmentFinal=(*it3).first;
            }
            it3++;
        }
        size_t found=(*it).first.find_last_of(".");
        string tTemplateName=RDPQuery.get_RDPID((*it).first.substr(found-8,8));      
        if(AssignmentFinal!=GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel))
        {
            Assignment<<count<<"\t"<<(*it).first<<"\t"<<GenBank.get_Taxonomy(tTemplateName, TaxonomicLevel)<<"\t"<<AssignmentFinal<<"\t"<<GenBank.get_Taxonomy(BestHits[(*it).first], TaxonomicLevel)<<"\n";
            map<string, float>::iterator it4=tAssignment.begin();
            while(it4!=tAssignment.end())
            {
                Assignment<<(*it4).first<<"\t"<<(*it4).second<<"\n";
                it4++;
            }
            count++;
        }                           
        tAssignment.clear();
        it++;
    }// while it
    cout<<"Assignment Test Done\n";
    return 1;
}
